function [Ypred] = lmm(data, label, k, s, lam)
%LMM function [Ypred] = lmm(data, label, k, s, lam)
%       larg-margin manifold classifier
%       uses fast maximum variance unfolding to
%       extract kernel
%
%       version 0.1, 04/11/2014
%       Suqi Liu

X = double(data);
Y = double(label);
n = size(X,1);
l = size(Y,1);

[IDX,D] = knnsearch(X,X,'K',k);

Ivec = reshape(repmat((1:n)',1,k-1),[],1);
IDXvec = reshape(IDX(:,2:k),[],1);
clear IDX;
Dvec = reshape(D(:,2:k).^2,[],1);
clear D;
Dvec = Dvec/max(Dvec);

A = sparse(Ivec,IDXvec,Dvec);
A = A + A';
A = (A~=0);

%check connected
[S,~] = graphconncomp(A,'Directed',false);
if(S > 1) 
    error('Graph not connected! Try more neighbors.');
end

B = adj2inc(A);
clear A;

L = double(B)'*double(B);
clear B;
DD = diag(L);
L = 2*sparse(1:n,1:n,DD) - L;
clear DD;

[V,~] = eigs(L,s+1,'SM');
clear L;
V = V(:,1:s);
Q = V(Ivec,:) - V(IDXvec,:);
clear Ivec;
clear IDXvec;

P = repmat(Y,1,s).*V(1:l,:);

cvx_begin sdp
cvx_solver sedumi
    variable t
    variable gam(s,1)
    variable M(s,s) symmetric
    minimize (t-lam*trace(M))
    subject to
        sum(Q*M.*Q,2) <= Dvec
        M > 0
        P*gam >= ones(l,1)
        [M gam; gam' t] >= 0
cvx_end

Ypred = int8(sign(V(l+1:n,:)*gam));